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ABSTRACT 

Gravitational Waves (GWs) are tiny ripples in the fabric of space-time predicted by Einstein's General Rel- 
ativity. Pulsar timing arrays (PTAs) are well poised to detect low frequency (10~ 9 - 10~ 7 Hz) GWs in the near 
future. There has been a significant amount of research into the detection of a stochastic background of GWs 
from supermassive black hole binaries (SMBHBs). Recent work has shown that single continuous sources 
standing out above the background may be detectable by PTAs operating at a sensitivity sufficient to detect 
the stochastic background. The most likely sources of continuous GWs in the pulsar timing frequency band 
are extremely massive and/or nearby SMBHBs. In this paper we present detection strategies including various 
forms of matched filtering and power spectral summing. We determine the efficacy and computational cost 
of such strategies. It is shown that it is computationally infeasible to use an optimal matched filter including 
the poorly constrained pulsar distances with a grid based method. We show that an Earth-term-matched filter 
constructed using only the correlated signal terms is both computationally viable and highly sensitive to GW 
signals. This technique is only a factor of two less sensitive than the computationally unrealizable optimal 
matched filter and a factor of two more sensitive than a power spectral summing technique. We further show 
that a pairwise matched filter, taking the pulsar distances into account is comparable to the optimal matched 
filter for the single template case and comparable to the Earth-term-matched filter for many search templates. 
Finally, using simulated data optimal quality, we place a theoretical minimum detectable strain amplitude of 
h > 2 x 10~ 15 from continuous GWs at frequencies on the order ~ 1 /T a ^ s . 
Subject headings: 



1. INTRODUCTION 

Low frequency (10~ 9 -10~ 7 Hz) GWs are expected from su- 
permassive black hole binary systems (SMBHBs), cosmic 
strings, the big bang and inflationary era of the early uni- 
verse. GWs from these sources can manifest themselves in 
different ways. Single nearby SMBHBs can produce resolv- 
able waves with periods on the order of years ( Wyithe & Loeb| 
[2003j|Sesana et al.|2009||S"esana & Vecc hio 2010). SMBHBs 
and cosmic stri ngs can also produce GW bur sts (Damour & 
|Vilenkin|[200T] |Siemens etaLl[2007i |Leblond et al.||2009T ln 
which the duration of the GW signal is much less than the 
observation time. We also expect PTAs to be sensitive to a 
stochastic background of unresolvable sources. Pulsar timing 
arrays (PTAs) offer an opportunity to detect low frequency 
GWs from all of these sources. The concept of a PTA com- 
posed of the best timed MSPs was firs t develo ped over two 
decades ago (Romani 1989; Foster & Backer 1990). Today 
there are three main PTAs in existence with the goal of GW 
detection using pulsars: the European Pulsar Timing Array 
(EPTA; [Janssen et al.|2008| l, the North American Nanohertz 
Observatory for Gravitational waves (NANOGrav; |Jenet et al. 
2009|l, and the Parkes Pulsar Timing Array (PPTA; Manch- 
ester 2008), all of which are in collaboration to form the In 



ternational Pulsar Timing Array (IPTA; pobbs et aL] |2010). 
There is also a large international effort for the construction 
of future generation radio telescope arrays such as the Square 
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Kilometer Array (SKA; Lazio 2009 ), for which a primary sci- 
ence goal will be GW astrophysics. 

Previous authors have developed statistical data-analysis 
methods fo r searches of the PTA data sets for s tochastic back- 
grounds (Jenet et al. 2005 1 Anholm et al. 2009; van Haasteren 
[etaT1|2009[ |Yardley et al.||2011[ ID emorest et al.||2012| )~and 
burst sources ^Finn & L ommen 2010J! However, studies 
into continuous GW detection have been more theoretical 
or "proof-of-principle" in nature, as opposed to a more rig- 
orous detection method aimed at real data-analysis pipeline 



implementation. Prior to the establishment of PTAs, Jenet 
|et al.| ( 2004) 1 used existing pulsar data to rule out the pro- 



posed SMBHB system 3C66B, a possible source of continu- 
ous GWs. This work looked for the signature of a continuous 
GW in real pulsar data through the use of Lomb-Scargle pe- 
riodograms and sug gested a metho d for directed searches of 
known sources. Yar dley et al.| ( |2010"| ) also relied on the Lomb- 
Scargle periodogram to determine the sensitivity of a PTA 
to continuous GW sources as a function of GW frequency. 
|Sesana & Vecchio] ( 20 1Q\ developed a Bayesian framework 
for the detection of continuous GWs from monochromatic 
SMBHBs in circular orbits. This work only included the earth 
term in the GW signal model and estimated the uncertainties 
one would expect on search parameters via the Fisher Infor- 
mation matrix, which is known to perform well in t he high 
signal-to-noise r atio (SNR) regime ( Vallisne ri|2008[ l. [Corbin 
|& Cornish] ( |2010| l have developed a Bayesian Markov Chain 
Monte-Carlo (MCMC) technique for parameter estimation of 
an evolving SMBHB system in which the pulsar term is taken 
into account in the detection scheme. This work took advan- 
tage of a signal model in which the GW frequency evolves 
signifi cantly in tim e to determine the pulsar distance. Most re- 
cently, |Lee|etaD(|20TT|l have developed parameter estimation 
techniques based on vector Ziv-Zakai bounds incorporating 
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the pulsar term and have placed limits on the detectable am- 
plitude of a continuous GW. In their work, a future PTA with 
SKA sensitivity is assumed, resulting in high SNR signals for 
their parameter estimation studies. All of these methods are 
promising for parameter estimation in a relatively high SNR 
(e.g. SNR=20 in |Corbin & Cornish| ( p010) l) limit and a favor- 
able signal modeL 

The aim of this paper it to determine the most sensitive 
practical detection technique for continuous GW sources in 
the PTA data making no assumptions about the SNR of the 
signal. This is done by comparing multiple detection tech- 
niques, using the minimum detectable amplitude as our figure 
of merit. The paper is ordered as follows. In section J2| we 
introduce the formalism and notation that we will use in the 



paper. In section 2.1 we define the GW signal from a SMBHB 
and der ive t he resulting GW induced pulsar timing residuals. 
Section |2T2] reviews methods of matched filtering, maximum 
likelihood detection techniques and power spectral summing. 
In section [3] we describe the simulated data sets that are used 
in this work. Section [4] describes the different detection tech- 
niques and discusses the main results of the paper. Finally, in 
section [5] we summarize our work and mention prospects for 
future work. 

2. METHOD 

2.1. PTA Response to a Continuous GW 

While PTAs are poised to detect a stochastic GW back- 
ground due to SMBHBs in the next five years, single, resolv- 
able sources ma y also be detec ted at expected five-year sensi- 
tivity limits ( [Sesana et al.|2009] >. A GW is defined as a metric 
perturbation to flat space time, 



habit , 0) = e + ab (Q)h + (t , SI) + e*(Sl)h x (t, SI), 



(1) 



where Si is the unit vector pointing from the GW source to 
the Solar System Barycenter (SSB) and h+, h x and e A h (A = 
+, x) are the polarization amplitudes and polarization tensors, 
respectively (See the Appendix for more details). The GW 
will cause a fractional shift in frequency, v, that can be defined 
by a redshift in the times-of-arrival (TOAs) 



5v(t,Sl) 



= -e A b (si)\- p -^Ah A (t 1 si), 
2 i+fip 



where 



Ah A (t,n) = h A (t e )-h A (t p ). 



(2) 



(3) 



Note that we use the standard Einstein summing convention. 
Here t e and t p are the time at which the GW passes the earth 
and pulsar respectively and p is the unit vector pointing from 
the SSB to the pulsar. Henceforth, we will drop the subscript 
"e" denoting the earth time unless otherwise noted. From ge- 
ometry we can write 6 



t p = t-D(l - cos fj,), 



(4) 



where Si - p = -cosfi and D is the distance to the pulsar. To 
simplify our notation we introduce the pulsar "antenna pattern 
functions" 



f a (si) = - 



1 



p"p" 



2 i+si-p 

6 Note that we use units in which c = G= 1 . 



(5) 



where the redshift is now written as 



5u(t,Sl) 



■■F A (ti)Ah A (t,n). 



(6) 



For this work we will only consider circular, non-precessing, 
monochromatic SM BHB systems, as the se are expected to be 
the most prominent (Sesa na et al.|2009] ). Astrophysica l justi- 
fication for the abo ve approximations can be found in Sesana 
|& Vecchio] ( |2010| l. The word monochromatic indicates that 
the frequency evolution of the system is slow enough that 
we can make the approximation that f(t p ) = f(t e ) = const. It 
should be noted that we use the observed redshifted values. 
For example, the chirp mass and frequency in the rest frame 
are M r = A4/(l+z) and f r = /o( 1 + z), respectively, where z 
is the cosmological redshift. Assuming a monochromatic sys- 
tem with a circular orbit and an orbital angular frequency of 
oj or b = 2(27r/o) = 47r/o, where /o is the GW frequency we can 
now write the polarization amplitudes as 



h+(t) = h [(1 +cos 2 t)cos20„cos(27r/of — 4>o) 
- 2 cos 1 sin 2(f>„ sin(27r/of — </>o)] 

h x (t) = h[(l +cos 2 t)sin20 n cos(27r/ o f-0o) 
+ 2 cos tcos2(/>„ sin(27r/of - </>o)] , 



(7) 



(8) 



where 4>q is the orbital phase of the binary at t = 0, t is the 
inclination angle, and cf)„ is the angle to the line of nodes. We 
define the amplitude h, 



h = 2 



M^(nf ) 2 ^ 



D, 



(9) 



The magnitude of the polarization amplitudes are proportional 
to the SMBHB chirp mass M = (MiM 2 ) 3/s /(Mi +M 2 ) 1/5 , the 
comoving distance to the source D c , and the GW frequency 
/o (twice the orbital frequency), and can be written as 



h ~ 8 x 10" 



./'« 



/ M 
\WMq 

2/3 



5/3 



Dc 
lOOMpc 



5x 10" 8 Hz 



(10) 



We now use the redshift given in Eq. [6] to compute the GW 
induced pulsar timing residuals 



rit^f^dt 



1 



F A (Cl)Ah A (t,(l) 



Wo 



F (SI) 



2tt/ 

h A (t e )-h A (t e -D(l+Sl-p)) 



(11) 



Here we have written out the explicit dependence on the pul- 
sar distance and sky location. The first term in square brack- 
ets is the so-called "earth term" because it refers to the metric 
perturbation at earth and is correlated in all sets of pulsar tim- 
ing residuals. The second term in the square brackets is the 
so-called "pulsar term" because it refers to the metric pertur- 
bation at the pulsar and is uncorrected in all sets of pulsar 
timing residuals. Notice that the pulsar term carries a depen- 
dence on the GW sky location SI, and all of the dependence on 
the pulsar distance, D. As with the amplitude, we can express 
the approximate amplitude of the GW induced timing resid- 
uals as a function of the SMBHB source parameters (Sesana] 
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etal. 2009), 



r ~ 25.7ns 



/ M \ 5/ V D. 



fo 



\io 9 MqJ V 100M p c 

1/3 



(12) 



s 5x 10" 8 Hz / 

The residuals can be written as a function of 8 +M parameters 

\ = {0,<t>,M,D c Jo,L,<t> n ,<i>o,D}, (13) 

where D is a vector of the M pulsar distances. Detecting and 
characterizing a signal that is a function of many parameters 
can be quite difficult and will be discussed in future papers. In 
this work we aim to give a baseline to the problem, in that, we 
will assume that all parameters are known and simply assess 
how well a particular detection method can confidently detect 
the signal. 

2.2. Matched Filtering and Power Spectral Summing 

Here we will outline matched filtering in terms of the log- 
likelihood, and a power spectral summing technique. First we 
will review matched filtering basics in the context of a PTA. 
The problem of d etecting a signal in noisy data is well studied 
(Wain stein & Zubakov|[T9 62). We assume that the noise in 
each pulsar is additive, stationary and gaussian. For this case, 
the data for each set of pulsar timing residuals Jc Q (f,) = x, a can 
be written as 

Xfa — ri a -\-n[ ai (14) 

where r la = r a (f,) and n la = n a (ti) are the signal and the noise 
in each data set. Here i refers to the time index and a refers to 
the pulsar number. As is the method in matched filtering, we 
want to compare our data to a signal template of known form. 
Here we define r ia = r a (tj, A) as our template of known form 
where A is the vector of search parameters given in Eq. 
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We define the inner product of two functions of time x(t,) 
and y(tt) as 

; ' (15) 



(x\y) = ^2x i (P) i jy j , 



hi 

where C is the covariance matrix of the noise. With this 
framework in place we can specialize to the case of our PTA 
data and templates with white noise. We find the inner prod- 
uct of our full set of residual data with the corresponding set 
of templates for that data set as 



(16) 



where er^ is the rms of the residuals from the ath pulsar. A 
Wiener optimal statistic can be defined as 



P(A): 



l*k(A)) 

(KA)|KA)) 



(17) 



If the noise is Gaussian and the signal is present, then the 
signal-to-noise ratio is given as 



(18) 



SNR=(p(A')) = V(KA')k(A')), 



where A' is the best estimate of the source parameters. For 
our data analysis purposes we use the log-likelihood as our 



matched filtering statistic. Under the assumption of gaussian 
noise, we define the likelihood function as the probability of 
the data x a {ti) given some set of model parameters A 



p(x|A) = C n0 rmeXp 



l -({x-r{\))\{x 



-r(A)) 



(19) 



where C norm is a normalization constant. We now define the 
relative likelihood as A(A) = p(x\X)/ p(x\Q), where p(x\Q) is 
the probability of the data given the null hypothesis. From 
this we define the log-likelihood function 



In A(A)= ( (*|r(A))-~(r(A)|r(A)) 



(20) 



By defining the log-likelihood in terms of the relative likeli- 
hood, we incorporate hypothesis testing (whether the GW sig- 
nal is present or not) and parameter estimation into one statis- 
tic. We can define the SNR in terms of the log-likelihood as 
follows 



(lnA(A')) = i(p(A')) 



1 



-SNR 



(21) 



To determine if a signal is present and to determine the 
source parameters A, one would need to search parameter 
space to find the maximum value of the likelihood. This 
can be done through grid based methods, Nested Sampling, 
or Markov Chain Monte-Carlo (MCMC). For this work we 
are only concerned with detection of a source. In this case, to 
claim a detection, our statistic (log-likelihood) must be greater 
than a threshold value determined by a specified false alarm 
value. 

While the log-likelihood has the ability to simultaneously 
carry out detection, parameter estimation and hypothesis test- 
ing, we now describe a method aimed at detection. In this 
method, we simply calculate the power spectrum of each set 
of pulsar timing residuals and then sum the power weighted 
by the variance of each data set and look for the maximum 
value over all frequency bins. We define our detection statis- 
tic as 

M S a (f) 



V = max 

a-l 



(22) 



where S a { f) and a 2 a are the one sided power spectrum and the 
variance of the ath pulsar data set, respectively. A detection 
is claimed when the value of V is greater than some threshold 
value Vo corresponding to a false alarm rate. 

3. SIMULATED PTA DATA SETS 

The simulated PTA data sets used for this work represent 
a best case scenario when it comes to data quality. While 
the quality of the data is optimal, the properties of the PTA 
(i.e. distances, sky location, rms, etc.) are meant to repre- 
sent a realizable case. The array consists of up to 100 pul- 
sars uniformly distributed in both the azimuthal angle <f> and 
in the cosine of the polar angle cos#. The pulsar distances 
are also drawn from a uniform random distribution ranging 
from 0.5-3 kpc. It should be noted that although there are 
known millisecond pulsars with distances less than 0.5 kpc, 
the simulated pulsar distances here are meant to represent an 
average for a typical PTA. The pulsar timing residuals are 
evenly spaced over 10 years with 250 TOAs for each pulsar 
simulating roughly bi-monthly observations, and rms values 
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drawn randomly from a uniform sample ranging from 100- 
300 ns. The noise is simulated to be white, gaussian, addi- 
tive, and stationary. In real pulsar timing data, the residu- 
als will be unevenly sampled and the noise may have various 
red components. Fortunately, recent work suggests that most 
NANOGrav appear to be mostly white with little to no red 
noise contributions (Perr odin et al.| 2012; Ell is et al.|2012| l In 
addition, the pulsar timing residuals will not be stationary, as 
a quadratic must be fit out of the data to account for the spin- 
down of the pulsar. Specifically, our definition of the inner 
product in Eq. [16] no longer holds as we will need to include 
the covariance matrix of the data and incorporate a linear op- 
erator that takes into account this fitting. These are issues that 
will need to be addressed in order to make a fully functional 
data analysis pipeline for continuous GW searches and will be 
addressed in future papers. However, here we will deal with 
the simple case to illustrate the efficacy of the studied search 
techniques on a data set of optimal quality. 

4. ANALYSIS 

In this section we compare four different detection tech- 
niques and determine their efficacy in terms of a minimum 
detectable amplitude as a function of the number of pulsars in 
the array. The four detection methods are the full matched fil- 
ter, the earth-term only matched filter, the pairwise matched 
filter, and a simple power spectral summing technique. For all 
of these methods we will assume that we know the parameters 
of the source exactly and will only be interested in the lowest 
amplitude that each method can detect. Though unrealistic in 
most astrophysical scenarios, this gives us a simple baseline 
for comparison of the four methods. 

4.1. Detection Methods 

The full matched filter includes both the earth and pulsar 
terms from Eq. fTTIand is thus a coherent search technique. 
This is the optimal detection statistic for a continuous wave 
buried in gaussian noise. However, this method has the ma- 
jor drawback that it is computationally expensive to carry out 
in practice as the pulsar distances must be added as search 
parameters. For example, if we have an array of 20 pulsars 
and want to search over just 100 trial distances for each pul- 
sar, then we will need to use at least 10 40 templates. If one 
does not use grid based methods, this number will drop sig- 
nificantly. However, even using more advanced methods like 
Nested Sampling or MCMC, including the pulsar distance as 
a search parameter is still computationally expensive. In the 
low SNR regime, where the likelihood surface is relatively 
flat and noisy, it may be impossible to get an accurate pulsar 
distance estimate with finite computational resources. 

The Earth-term-matched filter uses templates that only de- 
pend on the coherent earth term and treats the pulsar term as 
a noise source. The major advantage of using only the earth 
term from Eq. [Tllfor this detection method is that one does 
not need to include the pulsar distance in the search, making 
this method much less computationally expensive. Since we 
are not including the pulsar term in the analysis, we will al- 
ways measure an SNR that is lower than the injected value. 
Also, there is a strong correlation between the pulsar distance 
and the sky location of the GW source. This will cause the re- 
covered sky location to be biased and have larger uncertainties 
than the full matched filter case. This is illustrated in Figure 
[T]in which a large GW (SNR=1000) is injected into the data 
from 40 pulsars. Both the full matched filter and the Earth- 
term-matched filter search over sky location, and a sky map is 
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FIG. 1 . — Skymaps created using the full matched filter (top) and the Earth- 
term-matched filter (bottom). The injected signal was very large (SNR=1()00) 
for illustration purposes. The "+" symbol indicates the injected sky location. 
We can see that the sky localization is biased and has a large uncertainty for 
the Earth-term-matched filter as compared to the full matched filter. In this 
particular case we observe a 46% loss in recovered SNR (See online version 
for color figures). 

created where the color scale indicates the log-likelihood. We 
can see that the full matched filter does a good job of local- 
ization whereas the Earth-term-matched filter is biased with a 
much larger error box in the sky. For this realization, using 
the Earth-term-matched filter also results in a 46% decrease 
in SNR compared to the full matched filter. 

The pairwise matched filter is a method that will allow us to 
take advantage of the pulsar term without the hindrance of an 
unworkable number of templates. This is done by construct- 
ing the full matched filter for each pair of pulsars, including 
the distance as a search term, and then adding the likelihoods 
as defined in Equation |20]in a pairwise fashion 

In A pw = $>lKA))a0-^(*|KA))a0, (23) 

a</3 

where the a/3 subscript denotes an inner product using two 
pulsars and the sum indicates a sum over all M(M - l)/2 
unique pulsar pairs. While we still need the same number 
of templates for the intrinsic parameters of the source, we 
need far fewer distance templates. For example, in the case 
of M pulsars and 100 trial distances for each pulsar, the full 
matched filter requires at least 10 2M distance templates. How- 
ever, for the pairwise matched filter we only perform a co- 
herent search each pulsar pair, therefore for each pair we only 
require 10 4 distance templates making the total number of dis- 
tance templates M(M- l)/2 x 10 4 , which is significantly less 
then the full matched filter. This method is still negatively af- 
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fected by the strong degeneracy between pulsar distance and 
sky location but by using many pulsar pairs, this degeneracy 
will be greatly reduced. 

The most simple and computationally inexpensive method 
that one could use to detect a continuous GW is power spec- 
tral summing as described in Eq. |22] A method very similar to 
this was used in |Yardley et al.| ( |2010[ l on real data to produce 
sensitivity curves for the PPTA pulsars. This method is rel- 
atively robust in that it does not depend on any signal model 
templates. The disadvantages of this method are that it is in- 
coherent because it does not keep track of phase information 
and it gives no indication of the true parameters of the source. 
With real data that is irregularly sampled and may contain 
red noise processes, going into the Fourier domain may pose 
problems that will be addressed in future work. However, for 
this analysis we are using this method as a baseline to com- 
pare the matched filtering statistics. 

4.2. Efficacy of Detection Statistics 

Here we will outline our Monte-Carlo simulations and 
present our results. A good figure of merit for a detection 
statistic is the minimum detectable amplitude, that is, we wish 
to find the amplitude that can be detected above some thresh- 
old in 95% of the simulated realizations. To find the minimum 
detectable amplitude, we must first define a false alarm rate, 
that is, that rate at which we expect to make a detection when 
no signal is present. Assuming gaussian statistics, we want a 
4er detection significance corresponding to a false alarm rate 
of 1/15,787. All four of our detection methods will have a dif- 
ferent threshold value corresponding to this false alarm rate. 
To calculate these thresholds we perform the following simu- 
lation. We choose a template at random for each realization 
and calculate the detection statistic for 15,787 realizations of 
noise and record the maximum value. Statistically, this means 
that we would expect to get a value for our statistic that is 
above this maximum value ~ 0.016% of the time if our data 
was pure gaussian noise. 

A Monte-Carlo simulation was run to determine the mini- 
mum detectable amplitude as a function of the number of pul- 
sars in the array for four detection statistics. The steps of the 
simulation are as follows: (J) Choose the number of pulsars in 
the array and simulate residual data with white noise, (if) Fix 7 
the chirp mass A4, frequency /o, and distance D c to construct 
the amplitude h given in Equation|9] Then create a GW source 
in the sky with a given set of parameters X = (9,4>,t, 4> n ) drawn 
from random distributions. (Hi) Add the GW induced residu- 
als into the simulated PTA data, (iv) Run the detection statis- 
tic code assuming all parameters are known exactly 8 and out- 
put the log-likelihood. If the log-likelihood is above a given 
threshold value, count as a detection, otherwise, count as a 
non-detection, (v) draw a new parameter vector A and repeat 
steps Hi and iv. Repeat this for 10,000 GW source realizations 
(different realizations of A) and record the percent of sources 
detected, (vz) Keep M. and /o fixed and change D c to obtain 
a new amplitude h, repeat step v until 95% of the realizations 
are detected, (vz'z) Change the number of pulsars in the array 

7 For this work we fix the frequency to the lowest detectable frequency of 
fo = 1 /Tibs an d the chirp mass to a reasonable value of M = 5 X lO 8 A/0 

8 This step is somewhat different for the Earth-term-matched filter because 
the largest SNR does not correspond to the case where the filter signal param- 
eters are the same as the input signal parameters since we are not including 
the pulsar term. For this case we carry out a search over source sky location 
for each iteration in order to obtain the maximum possible log-likelihood. 



and repeat the entire procedure. In practice, a bisection root 
finding method is used to determine the detection probability 
instead of linearly increasing h until the 95% level is reached. 

4.2.1. Single Template Case 

We ran the simulation described above for each of our de- 
tection methods for PTAs with 15-100 pulsars. In this case 
we assume that there are no search templates except for the 
one exactly matching the data. This means that we will ob- 
tain the lowest possible false alarm probability. The results 
are shown in Figure [2] It is obvious that that the unrealizable 
full matched filter (thin dashed line) can detect the lowest am- 
plitudes of the four detection methods tested. However, if one 
uses the pairwise matched filter (thin dotted line), very little 
sensitivity is lost for PTAs with up to 30 pulsars. Also, by 
using only the earth term in a matched filter search (thin solid 
line), the resulting minimum detectable amplitude is only a 
factor of two higher than the optimal method of filtering for 
the entire signal. It is also important to note that these two 
matched filtering methods (full and earth term) scale roughly 
the same with the number of pulsars in the array, therefore, 
this factor of two is independent of the number of pulsars. 
The incoherent power spectral summing method does approx- 
imately four times worse than the optimal matched filter case 
at 20 pulsars in the array and will continue to get worse as the 
number of pulsars in the array increases since the trend has 
a weaker dependence on the number of pulsars. It is known 
that the SNR scales as \J~M for coherent methods and scales 
as (M) 1 ' 4 for incoherent methods, where M is the number of 
"detectors". The pairwise matched filter is also incoherent so 
the SNR has a flatter slope vs. N than when using the other 
two matched filtering methods. For a real PTA, the SNR will 
not scale exactly as mentioned above because each set of pul- 
sar timing residuals do not have the exact same characteris- 
tics (different noise, sky location, distance, etc.). As a sanity 
check, we ran the simulation on PTAs where every pulsar had 
the same sky location, distance and rms residual. In this case 
the curves in Figure [2] scale exactly as mentioned above. 

4.2.2. Multiple Template Case 

The above section deals with the case of one template. In 
reality, the act of searching over many templates will serve 
to increase the false alarm probability. In the case of gaus- 
sia n noise the false alarm rate can be calculated analytically 

as ( |Maggiore|2007| l 



/>fa : 



dpe 



-p-jla 



■■ 2erfc(/9o/v2er) 



(24) 



where po is a threshold value of the SNR, p is the output of 
the Wiener Filter in Eq. [17] a is the standard deviation of 
the probability distributionfunction and erfc(z) is the com- 
plementary error function. Formally, this false alarm rate is 
derived based on a Wiener filter defined in Eq. [T7j however, 
since the expectation value of the log-likelihooa is propor- 



tional the expectation value of the Wiener filter (see Eq. 21 
this equation is still valid for the log-likelihood method mat 
we use here. Because we deal with only white gaussian noise 
in this paper, our simulated thresholds are in agreement with 
Eq. [24] However when calculating a false alarm probability 
for a search with N templates, the total false alarm probability 
is Ppa = NppA, when /?fa is much less than unity. This implies 
that if one wants to keep the same detection significance (less 
than 1/15,787 chance of occurring in noisy data), then the 
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FIG. 2. — Plot of minimum detectable amplitude vs. number of pulsars for the full matched filter (thin dashed line), Earth-term-matched filter (thin solid line), 
power spectral summing (thick dash-dot line) and the pairwise matched filter (thin dotted line). Also plotted are the full matched filter (thick dashed line) , the 
Earth-term-matched filter (thick solid line), and the pairwise matched filter (thick dotted line) with a realistic number of search templates. To make these plots, a 
Monte-Carlo simulation was run to find the amplitude at which 95% of GW source realizations were detected for a given number of pulsars (See online version 
for color figures). 



threshold value must be increased. When this is done, we can 
see from Figure[2]that the minimum detectable amplitudes for 
the Earth term matched filter, pairwise matched filter and full 
matched filter are nearly the same because of the significantly 
smaller number of templates for the earth term and pairwise 
matched filters. It is also important to note the values of h on 
the y-axis. Since we are dealing with the best case scenario 
in terms of data quality (white gaussian noise, evenly spaced 
data, no timing model subtraction), this plot shows us that, at 
current levels of timing precision (rms ^100 ns), we could 
never confidently detect any signal with an amplitude below 
h = 9 x 10" 16 even with a PTA of 100 pulsars. At 20 pulsars, 
we could only possibly detect a source with h > 2 x 10~ 15 . 
In terms of placing limits on the minimum detectable ampli- 
tude for real data, our work suggests that if one implements an 
Earth-term-matched filtering technique a s opposed to a power 
spectra summing technique as used in Yar dley et al.| ( |2010| l, 
one could improve the results by a factor of ~ 2. 

5. SUMMARY 

In this work we have tested the efficacy of four detection 
techniques on simulated data sets when searching for con- 
tinuous GW signals from SMBHBs. We have shown that 
a matched filter using only the correlated earth term in the 
search templates results in a minimum detectable amplitude 
that is 2 times higher than the optimal matched filter using 
both the correlated and uncorrected terms. We have also 
shown that when using a pairwise matched filter it is possi- 
ble, in principle, to obtain nearly the same sensitivity as the 
full matched filter. However, when performing a real search, 
the strong correlations between sky location of the GW source 



and pulsar distance may cause problems with SNR recovery, 
a problem that does not affect the Earth-term-matched filter. 
We have also shown that an incoherent power spectra sum- 
ming method results in a minimum detectable amplitude that 
is 4 times higher at the present case of a 20 pulsar PTA and 
will continue to get worse as the number of pulsars in the ar- 
ray increases. When the number of search templates is taken 
into account, we find that using an Earth-term-matched fil- 
ter is nearly as sensitive as the full matched filter. Moreover, 
this work gives an idea of the prospects of detecting a con- 
tinuous GW with PTAs, by placing lower limits on the de- 
tectable amplitude for data of optimal quality. The advantages 
and disadvantages of the various detection methods have been 
discussed and it has been shown that using a full matched fil- 
ter with the pulsar distances included as search parameters is 
very computationally expensive (maybe even impossible for 
some cases). Because of this and the relatively low computa- 
tional cost along with increased sensitivity over power spec- 
trum techniques, the Earth-term-matched filter and pairwise 
matched filter appear to be practical choices for a detection 
method in a data analysis pipeline for use on real pulsar tim- 
ing data. 

This work gives some insight into what detection tech- 
niques should be used in a fully functional pipeline. More 
sophisticated data analysis methods will have to include the 
effects of irregularly sampled data, red noise (both correlated 
noise in the form of the stochastic GW background and un- 
correlated noise in the form of intrinsic timing noise and in- 
terstellar medium effects), and timing model parameter fits. 
The methods described here give basic detection algorithms 
that can be modified for use with real data. After a detection 
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is made, the next step is parameter estimation. This will re- in a large parameter space. Both of these challenges are cur- 
quire fast, efficient algorithms to find the correct parameters rently being studied and will be the subject of future papers. 

APPENDIX 

POLARIZATION TENSOR IN THE SSB REFERENCE FRAME 

Here we will show how one can convert the polarization tensors into the SSB coordinates. Once again, a GW is defined as a 
metric perturbation to flat space time, 

h ah (t,Cl) = e + ab (Cl)h + (t,Cl) + e^ h (n)h x (t,n), (Al) 

where f2 is the unit vector pointing from the GW source to the Solar System Barycenter (SSB) and h+, h x and e^ b {A = +, x) 
are the polarization amplitudes and polariza tion tensors, r espectively. The polarization tensors can be converted to the SSB by 
the following transformation. Following Wahlquist|( |1987| l we write the polarization tensors in terms of the wave principal axes 
described by unit vectors m and ft 

e + ab {Vt) = m a m b - h a h b , (A2) 
e^ b {Vt) = m a h b + n a m b . (A3) 

In the SSB coordinate center we define 

£1 = -(sinf?cos0)i-(sin6lsin0)3>-(cosf7)z, (A4) 
m = -(sin 0)i + (cos 4>)y, (A5) 
ft = -(cos 8 cos 0)i - (cos 8 sin (f>)y + (sin 8)t (A6) 

In this coordinate system, 8 = tt/2- S and = a are the polar and azimuthal angles of the source, respectively, where S and a are 
declination and right ascension in usual celestial coordinates. In this coordinate system, the polarization tensors can be written as 

sin 2 0- cos 2 cos 2 8 — cos 2 0(1+ cos 2 8) sin cos^cos^sint 
e + (8,4>)= | -cos 0(1 + cos 2 8) sin cos 2 0- cos 2 8 sin 2 cos 8 sin 8 sin | (A7) 
cos cos 8 sin 8 cos 8 sin sin 8 -sin 2 ( 

/ cos0sin20 -cos20cos( 
: (0,0) = -COS20COS6 1 -2 cos cos 8 sin cos0sinf? |. (A8) 
V — sin0sin0 cos sin L 
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